In [1]:
import pandas as pd
import numpy as np
import scipy.optimize 
import plotly.express as pl
In [161]:
import plotly.io as pio
pio.renderers.default='notebook'

Вопрос 1¶

Для первичного осмотра того, что происходит с retention, построим график retention в целом по всем клиентам.

In [162]:
data = pd.read_excel("test_data_1.xlsx")
data["order_date"]=data["order_date"].apply(pd.to_datetime)
data=data.sort_values(['order_date']).reset_index(drop=True)
In [163]:
data["time_diff"] = 0
dct = {i : 0 for i in data["client_id"].unique()}
In [164]:
for i in range(len(data)):
    if dct[data["client_id"].loc[i]] ==0:
        data["time_diff"].loc[i]=0
        dct[data["client_id"].loc[i]] = data["order_date"].loc[i]
    else:
        data["time_diff"].loc[i] = (data["order_date"].loc[i]-dct[data["client_id"].loc[i]])/np.timedelta64(1, 'D')
C:\Users\1\AppData\Local\Temp\ipykernel_13360\827410020.py:3: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

C:\Users\1\AppData\Local\Temp\ipykernel_13360\827410020.py:6: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

Я осознанно игнорирую первую точку на графике (по-сути, день первой покупки), потому что если её оставить график будет практически нечитаемый.

In [165]:
Y=(data[["time_diff","client_id"]].groupby("time_diff").nunique()/len(data["client_id"].unique())*100).iloc[1:]
fig = pl.line(Y)
fig.update_layout(
    xaxis_title="Days after first purchase", yaxis_title="Retention (%)"
)
fig.show("notebook")

Что мы видим? Во-первых, наблюдается существенная просадка через 30 дней после первой покупки, а также примерно через 11 дней. Если мы, скажем, подразумеваем, что речь идет про "Руки", то можно сказать, что такая ситуация согласуется со здравым смыслом по следующим соображениям:

  1. Большая часть ремонтных работ заканчивается уже на следующий день и в течение года может и вовсе не потребоваться дополнительных работ, а ко всему прочему, даже если работа занимает больше времени, чаще всего заказ требует только одного специалиста: только сантехника или только электрика. (Из собственного опыта: специалист из "рук" собирал мне кухонный гарнитру примерно в течение месяца, но транзакция то была одна!) - это обуславливает, что retention сразу падает до уровня менее 1%. Сюда можно отнести многие работы по обслуживанию сантехники, сборку одиночного элемента мебели, уставновку окна.
  2. Существует заметная, но не такая уж и существенная категория работ, укладывающаяся примерно в полторы недели: при работе над отдельной комнатой вполне может потребоваться несколько специалистов (и несколько транзакций): например, клиенту потребуется отдельно в сжатые сроки оплатить сборщика мебели, поклейщика обоев и электрика.
  3. Наконец, крайне редко в быту можно столкнуться с такой ремонтной работой, которая потребует работ более чем на месяц: это может быть, например, большой ремонт квартиры в целом. Само собой, в таком случаем мы столкнемся с несколькими подряд идущими транзакциями на разных специалистов, причем не исключено, что по несколько раз. Такие работы длительны, но редко могут затянуться больше чем на месяц.

Вооружившись этой информацией, проведем когортный анализ, где под когортой будем понимать клиентов, сделавших заказ в определенном месяце.

In [166]:
data['first_day'] = 0
for i in range(len(data)):
    data['first_day'].iloc[i] = dct[data["client_id"].iloc[i]]
C:\Users\1\AppData\Local\Temp\ipykernel_13360\2342964297.py:3: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

In [167]:
data['first_day'] = pd.to_datetime(data['first_day'], format='%m/%d/%Y %H:%M')
data['cohort'] = data['first_day'].dt.to_period('M')
In [168]:
data_cohort = data.groupby(['cohort', 'time_diff']).agg(n_customers=('client_id', 'nunique')).reset_index(drop=False)
In [169]:
pivot = data_cohort.pivot_table(index='cohort', columns='time_diff', values='n_customers')
In [170]:
cohort_size = pivot.iloc[:, 0]
retention_matrix = pivot.divide(cohort_size, axis=0).iloc[:, 1:]
In [171]:
fig = pl.imshow(retention_matrix.fillna(0).values,  aspect="auto",
                labels=dict(x="Days after first putchase", y="Cohort (month)", color="Retention (%)"),
                y=["January", "February"," March", "April", "May", "June", "July", "August", "September", "October"]
               )
fig.update_xaxes(side="top")
fig.show("notebook")

В целом, можем наблюдать аналогичную картину: после 30-го retention существенно проваливается. Что мы можем выяснить дополнительно?

Просадка на 10 (примерно на 10 я имею ввиду) день характерна не для каждого месяца: к ним относится январь, апрель, июнь, июль и август. Про апрель не слишком понятно, но остальные месяца - традиционные месяца отдыха, когда, вполне вероятно, люди не склонны планировать длительные и, что важнее, комплексные ремонтные работы.

Чем стоило бы дополнить анализ? Ну, судя по датасету, напрашивается посмотреть на более тонкие когорты: например, по номеру города. Ещё можно дополнительно сгруппировать клиентов по сумме заказа (дешевые, средние и дорогие).

В любом случае, напрашивается смотреть не только на удержание, но и на сумму заказов: возможно, что те крохи клиентов, которые остаюстя после месяца, могут приносить существенную выручку.

Вопрос 3¶

Решение обобщенным методом моментов¶

$$z=(a+x)^{b+y}+\epsilon$$

Для оценивания двух параметров необходимо, по меньшей мере, 2 моментных условия. Cоставим моментные условия следующим образом: $$M(z)=M((a+x)^{b+y})+0.05$$ $$M(z^2)=M((a+x)^{2b+2y})+0.1M((a+x)^{b+y})+\frac{0.01}{3}$$ Заменим моменты на выборочные аналоги и приведем моментные условия к виду $g_i=0$: $$g_1=\frac{\sum_{i=1}^{n} z_i}{n}-\frac{\sum_{i=1}^{n} (a+x_i)^{b+y_i}}{n}-0.05$$ $$g_2=\frac{\sum_{i=1}^{n} z_i^2}{n}-\frac{\sum_{i=1}^{n} (a+x_i)^{2b+2y_i}}{n}-0.1\frac{\sum_{i=1}^{n} (a+x_i)^{b+y_i}}{n}-\frac{0.01}{3}$$ Будем оценивать параметры, решая следующую задачу: $$g_1^2+g_2^2\rightarrow min$$

In [119]:
#проверка на синтетических данных
x = np.random.uniform(low=1, high=10, size=100)
y = np.random.normal(loc=1, scale=10, size=100)
eps = np.random.uniform(low=0.0, high=0.1, size=100)
a = 1
b = 9
z = np.float_power((a+x),(b+y))+eps
df = pd.DataFrame({"x":x,"y":y,"z":z})
#df
def estimator(x):
    a,b=x
    global df
    g1=df["z"].mean()-((a+df["x"])**(b+df["y"])).mean()-0.05
    g2=(df["z"]**2).mean()-((a+df["x"])**(2*(b+df["y"]))).mean()-0.1*((a+df["x"])**(b+df["y"])).mean()-0.01/3
    return g1**2+g2**2
rranges = (slice(1, 10, 0.01), slice(1, 10, 0.01))
x1 = scipy.optimize.brute(estimator,rranges)
x1
Out[119]:
array([1., 9.])
In [118]:
df_3 = pd.read_excel("test_data_3.xlsx")
def estimator(x):
    a,b=x
    global df_3
    df = df_3
    g1=df["z"].mean()-((a+df["x"])**(b+df["y"])).mean()-0.05
    g2=(df["z"]**2).mean()-((a+df["x"])**(2*(b+df["y"]))).mean()-0.1*((a+df["x"])**(b+df["y"])).mean()-0.01/3
    return g1**2+g2**2
rranges = (slice(1, 10, 0.01), slice(1, 10, 0.01))
x1 = scipy.optimize.brute(estimator,rranges)
x1
Out[118]:
array([3.34305807, 2.58991463])

Вопрос 2¶

In [235]:
df_2 = pd.read_excel("test_data_2.xlsx")
In [236]:
from sklearn.cluster import DBSCAN
In [237]:
clstr = DBSCAN(eps=0.1,min_samples=5,p=1)
clstr.fit(df_2)
clstr.labels_
Out[237]:
array([-1,  0,  1, -1, -1,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,  2,
       -1, -1,  3, -1, -1, -1, -1,  0,  3,  2,  1, -1, -1, -1, -1, -1, -1,
       -1, -1,  0,  4,  3, -1, -1, -1, -1,  5,  5,  2, -1,  0, -1, -1,  1,
        4, -1, -1, -1,  3, -1, -1, -1,  6, -1,  2, -1, -1,  5, -1,  6, -1,
        6, -1, -1, -1,  4, -1, -1, -1, -1, -1, -1,  5, -1, -1, -1,  4, -1,
       -1, -1,  0, -1, -1,  6,  6, -1, -1,  5, -1, -1,  4, -1, -1, -1, -1,
       -1,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,  1, -1,
        2,  7,  7, -1,  8, -1, -1,  9, -1, -1, -1, -1, -1, -1,  9, -1, -1,
       -1, -1, -1, -1, -1, -1, -1, -1,  9, -1, -1, 10, -1,  7,  7, -1, -1,
       -1, -1,  8, -1,  8, 11, -1, -1, 11, 10, -1,  9,  9, 10, -1, -1,  9,
        9, -1, 11, -1, -1, -1, -1, -1, -1, -1, -1, 10,  8,  8,  8, -1, -1,
       -1, -1, 11,  7, -1,  8, -1, -1, -1, -1, -1, -1, -1,  8, -1, -1, 10,
       11, -1, -1, -1, -1, -1, -1, -1, -1, 10, -1, -1, -1, -1,  9, -1, -1,
        8, -1, -1, -1, -1, -1, -1, -1, -1, -1,  7, -1, -1, -1, -1, -1, -1,
       -1, -1], dtype=int64)
In [238]:
#df_2["x"] = X
df_2["lbls"] = clstr.labels_
df_2.groupby("lbls").mean()
Out[238]:
0 1 2 3 4 5 6
lbls
-1 0.502924 0.415205 0.497076 0.538012 0.520468 0.526316 0.532164
0 0.000000 1.000000 0.000000 0.000000 1.000000 0.000000 1.000000
1 0.000000 1.000000 0.000000 0.000000 1.000000 0.000000 0.000000
2 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000
3 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000 1.000000
4 0.000000 1.000000 0.000000 1.000000 1.000000 0.000000 1.000000
5 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000 0.000000
6 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000
7 1.000000 1.000000 0.000000 1.000000 0.000000 0.000000 0.000000
8 1.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000
9 1.000000 1.000000 1.000000 0.000000 0.000000 0.000000 0.000000
10 1.000000 1.000000 1.000000 0.000000 1.000000 0.000000 0.000000
11 1.000000 0.000000 1.000000 1.000000 1.000000 1.000000 0.000000

Понятно, что алгоритм кластеризации склоняется к тому, чтобы выделить вершины гиперкуба (либо просто весь гиперкуб при другом параметре $\epsilon$), что, потенциально, не совсем верно (совсем не верно) т.к. чисто технически это скорее просто одинаковые точки, а не кластеры, то есть не подмножества "похожих" точек. К тому же, даже если придерживаться идеи кластеризации по принципу "кластер==вершина гиперкуба" при большей "выборке" можно вообще дойти до того, что кластеров там будет $2^7$, что не круто с любых возможных прикладных позиций и, вероятно, абсолютно теряет всякую интепретируемость.

Здесь я бы подумал о том, чтобы бинарные вектора привести к десятичным числам и кластеризовать уже их:

In [239]:
def binatointeger(binary):
    number = 0
    for b in binary:
        number = (2 * number) + b
    return number#
X=[binatointeger(df_2.loc[i].values) for i in range(len(df_2))]
In [258]:
from sklearn.metrics import silhouette_score
In [266]:
for i in [1,5,10,15,20,25]:
    clstr = DBSCAN(eps=i,min_samples=5,p=1)
    clstr.fit(np.array(X).reshape(-1, 1))
    lbl = clstr.labels_
    print(f" Силуэт: {silhouette_score(df_2[[1,2,3,4,5,6]],lbl)} eps: {i} n: {set(lbl)}")
 Силуэт: -0.049691614250458205 eps: 1 n: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, -1}
 Силуэт: 0.013905347897382419 eps: 5 n: {0, 1, 2}
 Силуэт: 0.033833562157288305 eps: 10 n: {0, 1}
 Силуэт: 0.033833562157288305 eps: 15 n: {0, 1}
 Силуэт: 0.033833562157288305 eps: 20 n: {0, 1}
 Силуэт: 0.033833562157288305 eps: 25 n: {0, 1}

Ну вот из соображений максимизации силуэта я бы предположил, что тут два кластера, которые предсттавляют из себя что-т вроде двух противоположных вершин гиперкуба, вокруг которых разбросаны точки.

Вопрос 4¶

Решение¶

Идея заключается в следующем: понятно, что диверсификация портфеля - это хорошее решение, поскольку позволяет найти разумный компромисс между доходностью и риском. В связи с этими соображениями, я бы разделил все доступные источники на 3 группы: гарантированные, консервативные и рискованные.

К группе гарантированных отнесем (индекс d), например, вклад по депозиту (теоретики обычно считают, что это ОФЗ). Этот актив характеризуется отсутствием риска и низкой дохожностью.

К группе консервативных (индекс g) отнесем активы, которые могут как падать в цене, так и сокращаться, то есть риск по ним имеется, но небольшой, зато они дадут большую доходность, чем предыдущая группа. Эту группу будем условно считать вложением в золото.

Наконец, рискованным активом (индекс s) будет вложение в акции. Для упрощения будем считать, что мы можем вкладываться только в фонды, например в S&P500. Понятно, что это наиболее рискованно, но и наиболее доходно.

Мы будем работать со средней доходностью портфеля, которую можно записать так:

$$R=w_dr_d+w_gr_g+w_sr_s$$

Где $r_i$ - средняя доходность i-го актива, а $w_i$ - доля i-го актива в портфеле

Соответсвенно, дисперсия (мера риска) доходности портфеля может быть записана так:

$$D=w_d^2\sigma_d^2+w_g^2\sigma_g^2+w_s^2\sigma_s^2=w_g^2\sigma_g^2+w_s^2\sigma_s^2$$

(здесь мы, для простоты, считаем, что активы нескоррелированы и доходность вклада по депозиту безрисковая, то есть случайной величиной не является)

Положим, взяв дисперсии с потолка, а доходности из даннных за последний год, что:

$$r_d=16 ;\sigma_d^2=0$$

$$r_g=20 ;\sigma_g^2=1$$ $$r_s=25 ;\sigma_s^2=2$$

Я консервативен, поэтому решать буду следующую задачу:

$$1w_g^2+2w_s^2 \rightarrow min$$

S.T. $$16w_d+20w_g+25w_s=22$$ $$w_d+w_g+w_s=1$$ То есть я хочу минимизировать риск по портфелю, соглашаясь при это на доходность в 22% годовых.

In [267]:
def disp(x):
    w1,w2,w3=x
    return w2**2+2*w3**2
In [277]:
x1=scipy.optimize.minimize(disp, x0=[1/3,1/3,1/3],
         constraints = ({'type': 'eq', 'fun': lambda x:  x[0] + x[1] + x[2]-1},
        {'type': 'eq', 'fun': lambda x: 16*x[0] + 20 * x[1] + 25 * x[2]-22}))
print(f"Итого я размещу {x1.x[0]*1000000} руб. на депозите, {x1.x[1]*1000000} руб. в виде золота и {x1.x[-1]*1000000} руб. вложу в акции")
Итого я размещу 97346.21039558554 руб. на депозите, 424776.8212879459 руб. в виде золота и 477876.9683164685 руб. вложу в акции

Вопрос 5¶

В моем понимании, ответ на этот вопрос в существенной степени сопряжен с конкретными задачами моделирвания. Скажем, если стоит задача полуения прозрачно интерпретируемой модели (например, линейной), тогда со статистической точки зрения нам нужна хорошая стратегия идентификации параметров и, в общем-то, хороший дизайн эксперимента. Скажем, в линейной модели пропуск так называемых существенных переменных приводит к несостоятельности оценок, что для интерпретируемой и, по-сути, описательной модели, крайне нежелательно. В этом смысле соц-дем данные нам скорее необходимы, чем просто "желательны", в то время как все прочие признаки мы могли бы попытаться использовать как, например, инструментальные переменные.

Если же решаемая задача носит скорее инженерный характер и от модели в первую очередь требуется давать адекватный прогноз на отложенной выборке, я бы подумал о том, как сопоставить возможные потери от неиспользования этих данных, поскольку их сбор, что очевино, тоже требудет каких-то ресурсов. Скажем, если данные добываются посредством опроса, следовало бы составить выборку (репрезентативную без этих признаков составлять сомнительно, наверное всё таки случайную), провести опрос, получить данные и сравнить метрики на отложенных выборках. Вероятно здесь также будет уместен подход feature importance, SHAP или какой-то другой, позволяющий оценить важность признаков для нашей модели: в конечном итоге может оказаться, что на самом деле самый важный признак - это какой-нибудь цвет кожи или этническая принадлежность, который можно и без привлечения доп ресурсов достать на основании информации об имени объекта.

В любом случае, моя исследовательская практика говорит о том, что соц дем данные могут быть существенно зашумлены и мы можем наблюдать неадекватные ответы, которые даны чисто формально в связи с нежеланием человека давать на них ответы. На практике это может приводить как к не очень страшным последствиям в виде урезания выборки, так и к очень серьехным, связанным с self selection. Последнее, вообще говоря, может существенно влиять на резульаты в контексте обоих подходов к моделированию.